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ABSTRACT 


We report the detection of pulsed gamma-ray emission from the fast millisecond 
pulsars (MSPs) B1937+21 (also known as J1939+2134) and B1957+20 (J1959+2048) 
using 18 months of survey data recorded by the Fermi Large Area Telescope (LAT) and 
timing solutions based on radio observations conducted at the Westerbork and Nangay 
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radio telescopes. In addition, we analyzed archival RXTE and XMM-Newton X-ray 
data for the two MSPs, confirming the X-ray emission properties of PSR B1937+21 
and finding evidence (~ 4cr) for pulsed emission from PSR B1957+20 for the first time. 
In both cases the gamma-ray emission profile is characterized by two peaks separated by 
half a rotation and are in close alignment with components observed in radio and X-rays. 
These two pulsars join PSRs J0034— 0534 and J2214+3000 to form an emerging class of 
ganuna-ray MSPs with phase-aligned peaks in different energy bands. The modeling of 
the radio and gamma-ray emission profiles suggests co-located emission regions in the 
outer magnetosphere. 

Subject headings: gamma rays: observation - pulsars: general - pulsars: individual 
(PSR B1937+21, PSR B1957+20) - radiation mechanisms: non-thermal 


Introduction 


The Large Area Telescope (LAT) aboard the Fermi Gamma-ray Space Telescope has firmly 
established millisecond pulsars (MSPs), rapidly-rotating neutron stars (P < 30 ms) with small 
rotational spin-downs (P < 10 -17 ), as sources of GeV gamma rays. Nine MSPs known prior to the 
Fermi mission have so far been observed to emit pulsed gamma rays (Abdo et al. 20 09d| a, 2010b), 
and radio searches at the position of Fermi LAT unassociated sources, such as those in the Fermi 


LAT First Year Catalog (Abdo et al. 

2010c 

), have led to the discovery of over thirty previously 

unknown MSPs (see e.g. 

Hessels et al. 

2011 

). Six of these new pulsars have been shown to emit 


has also detected gamma-ray emission from several globular clusters, and the observed properties 


are consistent with the summed contribution of a population of MSPs (Abdo et al. 2009c Kong 


et al. 2010 Abdo et al. 2010a). In addition, the AGILE telescope reported a 4.2cr detection of 
PSR B1821— 24 in the globular cluster M28 in gamma rays (Pellizzoni et al. 2009 ) . These different 
observations indicate that MSPs are prominent sources of gamma rays and that many of them are 
awaiting detection with the Fermi LAT. 

All MSPs detected by Fermi to date are relatively energetic, with spin-down luminosities 


E = Ait 2 IP/P 3 > 10 33 erg s 1 


(where I denotes the moment of inertia, assumed to be 10 45 g cm 2 


in this work), making PSRs B1937+21 ( E = 1.1 x 10 36 erg s 1 


and B1957+20 (E = 7.5 x 10 34 
erg s _1 ) good candidates for detection in gamma rays with Fermi. Nevertheless, the two MSPs are 
more distant than the bulk of gamma-ray-detected MSPs (see Abdo et al. 2010e), and are located 
at low Galactic latitudes and therefore suffer from strong contamination from the diffuse Galactic 
emission (the properties of these two MSPs are listed in Table [I]), making these pulsars difficult to 
detect. 


In this article we describe the detection of pulsed gamma-ray emission from PSRs B1937+21 
and B1957+20 using the first 18 months of data recorded by the Fermi LAT. In addition, we 
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analyzed the X-ray properties of the two MSPs using archival RXTE and XMM-Newton data. The 
two MSPs are observed to emit radio, X-rays, and gamma rays, in near alignment. We present 
results of the modeling of these radio and ganmra-ray components in the context of geometrical 
models of emission from pulsar magnetospheres. Additionally, we examine the possibility that 
gamma rays are produced by colliding winds in the PSR B1957+20 system, as was observed in 
X-rays (Stappers et al. 2003). 


2. PSRs B1937+21 and B1957+20 


PSR B1937+21 was the first MSP ever discovered (Backer et al. 1982), and remained the 
pulsar with the shortest known rotational period ( P ~ 1.558 ms) until the recent detection of a 
1.396 ms pulsar in the globular cluster Terzan 5, PSR J1748— 2446ad (Hessels et al. 2006). With a 


pulse period of 1.607 ms, the first ever “black widow” pulsar discovered, PSR B1957+20 (Fruchter 


et al. 1988), has the third shortest rotational period of currently known pulsars. The distances 


to these MSPs are relatively uncertain: the NE2001 model of Galactic electron density (Cordes 


& Lazio 2002) places PSR B1937+21 at d = 3.6 ± 1.4 kpc, assuming a 40% uncertainty in the 


dispersion model (Brisken et al. 2002). However, timing measurements by Verbiest et al. (2009) 
led to 7r = 0.13 ± 0.07 mas, and therefore d = 7.7 ± 3.8 kpc. The only distance estimates currently 
available for PSR B1957+20 are based on the dispersion measure (DM) and the models of Galactic 
free electron density. The NE2001 model places this MSP at d = 2.5 ±1.0 kpc. This estimate 
is supported by spectral analyses of the binary companion at optical wavelengths, placing a lower 


limit on the distance of d > 2 kpc (van Kerkwijk et al. 2011). In the following, we will use the 


parallax distance from Verbiest et al. (2009) of 7.7 kpc for PSR B1937±21 and the NE2001 distance 


of 2.5 kpc for PSR B1957±20. 

The proper motion of PSR B1937±21 is also relatively uncertain: available estimates b ased 
on pulsar timing measurements range from /r t = \J Ra cos 2 8 ± = 0.28 ±0.09 mas yr -1 (Verbiest 

2009) to [it = 1-6 ± 0.2 mas yr -1 (Hotan et al. 2006), while Kaspi et al. (1994), Cognard et al. 


(1995), and Verbiest et al. (2009) give intermediate values with comparably small uncertainties, 
possibly underestimated. However, even with the largest of the measurements, the period derivative 
is weakly affected by the kinematic Shklovskii effect, which makes the apparent P greater than 


the intrinsic value by {P^d)/c (Shklovskii 1970). At a distance of 7.7 kpc and assuming the 


largest proper motion value of 1.6 mas yr -1 , the expected Shklovskii contribution is 7 X 10 23 , 
negligible compared to the apparent P of ~ 1.05 x 10~ 19 . In this paper we will use the intermediate 
proper motion value measured by Cognard et al. (1995) of 0.80 ± 0.02 mas yr -1 . In contrast 


with PSR B1937±21, the apparent period derivative of PSR B1957±20 is strongly affected by the 
Shklovskii effect: with its distance of 2.5 kpc and transverse velocity of fix = 30.4 ± 0.6 mas yr -1 


(Arzoumanian et al. 1994), the Shklovskii effect decreases the P value by more than half. Table 
[l] lists characteristic properties of PSRs B1937±21 and B1957±20. The spin-down power E and 
magnetic held at the light cylinder Blc values are the highest among known Galactic MSPs. 
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These two pulsars have been extensively studied at high energies: in X-rays, PSR B1937+21 
has a two-peaked profile similar to its radio profile, with peaks in close alignment with the giant 


radio pulse emission regions, lagging the normal radio emission peaks slightly (Cusumano et al. 


2003). The X-ray emission from this pulsar is non-thermal, which distinguishes it from many other 


MSPs. The high E of PSR B1957+20 would seemingly make this pulsar a good candidate for 
detection of X-ray pulsations. However, the X-ray emission from the system clearly has a strong 
contribution from the interaction of the pulsar wind with the companion, which is modulated at 
the orbital period (Stappers et al. 2003). Previous searches for X-ray pulsations (see e.g. Huang & 


Becker 2007) have been unsuccessful. 


3. Observations and data analysis 


3.1. Radio timing observations 

The timing solutions predicting rotational and orbital behaviors of PSRs B1937+21 and 
B1957+20 as a function of time were obtained from radio timing measurements contemporane- 
ous with the first 18 months of the Fermi mission. With spin-down energies above 10 34 erg s _1 , the 
two MSPs have been monitored by radio telescopes around the world as part of the pulsar timing 
campaign for Fermi (Smith et al. 2008). For PSR B1937+21, we have built a timing solution by 


using 80 Times of Arrival (TOAs) taken at the Nangay radio telescope in France (Cognard et al. 


2011) at 1.4 GHz between 2008 May 27 and 2010 February 7 and with a mean uncertainty on the 


determination of individual TOAs of 45 ns, and 42 TOAs recorded at the Westerbork Synthesis 
Radio Telescope (WSRT) ( Voute et al.|2002 Karuppusamy et al.|2008 ) at 1.4 and 2.3 GHz between 
2008 January 27 and 2009 November 25, with a mean uncertainty of 579 ns. For PSR B1957+20, 
the timing solution was built using 38 TOAs recorded between 2006 May 14 and 2009 December 19 
at the Nangay radio telescope at 1.4 GHz with a mean uncertainty of 2.3 [as, and 426 WSRT TOAs 
recorded at 0.35 GHz between 2008 January 27 and 2010 January 30 with a mean uncertainty of 
2.8 /xs. 

Ephemerides were built using the Tempo2 pulsar timing packag^j] (Hobbs et al. 2006). The 
resulting ephemeris for PSR B1937+21 gives a Root Mean Square (RMS) of the timing residuals of 
197 ns. The Dispersion Measure (DM), necessary for the relative phasing of observations at different 
wavelengths, was fitted using the multi-frequency radio TOAs in the case of PSR B1937+21. We 
measured a DM value for this pulsar of 71.01931 ± 0.00019 pc cm -3 across the observation, where 


the error bar is the nominal la uncertainty reported by Tempo2. Cognard et al. (1995) measured 
a DM time derivative of —0.0012 ± 0.0001 pc cm” 3 yr . In our analysis we found no indication of 
significant time variation. In the case of PSR B1957+20, however, significant DM variations with 
time were observed and needed to be taken into account when building the timing solution. The 


1 http:/ /tempo2. sourceforge.net/ 
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dispersion measure for PSR B1957+20 was measured by generating one WSRT TOA per band of 
20 MHz for each observation, thereby producing 8 TOAs per observation. TOAs corresponding 
to weak detections or where the pulsar was in eclipse were discarded. Assuming that for a given 
epoch pulses in the different frequency sub-bands were emitted simultaneously, we fitted the DM 
by measuring the time drift of TOAs as a function of radio frequency. DM excursions of 0.0008 and 
—0.0011 around an average value of 29.1259 pc cm -3 were observed. A DM of 29.12644 ± 0.00029 
pc cm -3 was measured at the epoch TZRMJD of ~54550 MJD, the TZRMJD parameter defining 
a reference epoch at which the rotational phase predicted by the timing solution is 00 DM values 
for each TOA in our dataset were determined by interpolating the values measured with the WSRT 
data. The TOAs and their corresponding DM values were then analyzed using Tempo2, resulting 
in an RMS of timing residuals of 4.1 /rs. The latter RMS value is larger than the timing accuracy of 
the LAT of less than 1 /is (jAbdo et al. 2009e), but negligible for the low-statistics gamma-ray light 
curves of PSR B1957+20 (see Section 3.2). Similarly, the uncertainties on measured DM values 
of PSRs B1937+21 and B1957+20 correspond to errors in the extrapolation of 1.4 GHz TOAs to 
infinite frequency of 400 and 600 ns respectively, again negligible considering the low statistics. The 
timing solutions will be made available through the Fermi Science Support Centei[^} 


3.2. Gamma-ray analysis 


3.2.1. Initial searches for pulsations and spectral analysis 

To search for gamma-ray pulsations from the two MSPs, we selected events recorded between 
2008 August 4 and 2010 February 13, with energies above 0.1 GeV, zenith angles < 105° and 
belonging to the “Diffuse” class of events under the P6_V3 instrument response functions (IRFs), 


those events having the highest probabilities of being photons (Atwood et al. 2009). We excluded 
times when the rocking angle of the instrument exceeded 52°, required that the DATA.QUAL 
and LAT.CONFIG are equal to 1 and that the Earth’s limb did not infringe upon the Region of 
Interest (ROI). The gamma-ray data were analyzed using the Fermi science tools (STs) v9rl8p(0 
Photon dates were phase-folded using the Fermi plug-in distributed with the Tempo2 pulsar 
timing packag^] (Ray et al. 2011) and the ephemerides described in Section 3.1 For ROIs of 0.8° 
radii around PSRs B1937+21 and B1957+20, we obtained values of the bin-independent H - test 


parameter (de Jager & Biisching 2010) of 30.8 and 41.1 respectively, corresponding to pulsation 


significances of 4.6 and 5.4a. These relatively low significances, for PSR B1937+21 in particular, 
prompted us to verify the pulsed nature of the observed gamma-ray signal with an alternative 


2 See http://www.atnf.csiro.au/research/pulsar/ppta/tempo2/manual.pdf for more details 
2 http:/ /fermi. gsfc.nasa.gov/ssc/data/access/lat/ephems/ 

4 http:/ /fermi. gsfc.nasa.gov/ssc/data/analysis/scitools/overview.html 
5 See http://fermi.gsfc.nasa.gov/ssc/data/analysis/user/Fermi_plug_doc.pdf 
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As discussed in Kerr (2011), the LAT sensitivity to gamma-ray pulsars can be improved by 
weighting each photon by its probability of originating from the considered pulsar, and by taking 
these weights into account in the calculation of the pulsation significance. Among other advantages, 
this method is efficient at discriminating background events, which is particularly important for 
these two MSPs located at low Galactic latitudes and therefore observed in the presence of intense 
background contamination. In order to calculate the photon probabilities, we analyzed the spectral 
properties of putative gamma-ray sources located at the positions of the two MSPs. The spectral 
analysis was done using a binned maximum likelihood method ((Mattox et al. 1996) as implemented 
in the pyLikelihood python module in the Fermi STs. All point sources from the 1FGL catalog 
(Abdo et al. 2010c) within 15° of the radio location of PSR B1937+21 were included in the model, 
as well as additional sources from a list internal to the LAT team, based on 18 months of data. All 
point sources were modeled with power-law spectra except for the two known gamma-ray pulsars 
in the ROI, PSRs J1954+2836 and J1958+2846 (Abdo et al. 2009b; Saz Parkinson et al. 2010), 
which were modeled as exponentially cut off power laws, of the form: 


dN 

dE 


-r 


N 0 


(idv) exp 



(a 


In Equation (JTJ) , Nq denotes a normalization factor, V is the photon index, E c is the cutoff 
energy of the pulsar spectrum, while (3 is a parameter determining the steepness of the exponential 
cutoff. Spectra measured so far by the Fermi LAT for gamma-ray pulsars are generally well- 
described by simple exponential models, /3 = 1. The Galactic diffuse emission was modeled using the 
glLiem-v02 mapcube. The extragalactic diffuse background and residual instrument background 
were modeled jointly using the isotropicAem-v02 template]^ The normalizations and indices for 
all point sources within 8° of PSR B1937+21 were left free in the fit as well as cutoff energies E c 
for pulsars. The normalizations of the diffuse components were also left free. In a first attempt 
at analyzing the spectra of the two MSPs we considered all photons with energies above 0.1 GeV. 
However, PSR B1937+21 could not be detected with sufficient significance below 0.5 GeV, and for 
both pulsars the best-fit spectral parameters were strongly affected by the background emission. 
We therefore subsequently rejected the events with energies below 0.5 GeV. After a first iteration of 
the analysis, the 1FGL sources J1938.2+2125c and J1959. 6+2047, located 21.7’ and 0.5’ away from 
PSRs B1937+21 and B1957+20 respectively, were no longer found to be significant gannna-ray 
sources. We thus removed them from the spectral model and made another iteration. This led us 
to the gamma-ray spectral parameters listed in Table [2j where the first errors quoted are statistical 
and the second ones are systematic, and were calculated by following the same procedure as above, 


°Both diffuse models are available through the Fermi Science Support Center (FSSC) (see 
http: / /fermi. gsfc.nasa.gov/ssc/ 



but using bracketing IRFs for which the effective area has been perturbed by ±10% at 0.1 GeV, 
±5% near 0.5 GeV, and ±20% at 10 GeV with linear extrapolations in log space between. Table 
[2] also lists integrated photon fluxes F and energy fluxes G above 0.5 GeV. To enable comparison 
with previously observed gamma-ray MSPs we quote photon and energy fluxes extrapolated to 
lower energies, obtained by integrating Equation 0 for energies above 0.1 GeV. 

The gamma-ray spectra for PSRs B1937±21 and B1957±20 are shown in Figures [l] and [2j 
From Figure [2] it can be seen that the full energy range fit for PSR B1957±20 agrees well with 
the individual energy band fits. For PSR B1937±21 the full energy range fit is observed to be 
systematically below energy band fits between 0.5 and 1 GeV, which may suggest that the actual 
spectrum is softer than the one fitted using the entire energy range. Maximum likelihood fits over 
the entire energy range were also performed with the two MSPs modeled with simple power-law 
spectra (/3 = 0) in order to address the significance of measured cutoff energies. Using a likelihood 
ratio test we find that exponentially cut off power-law models (f3 = 1) are preferred at the 2.8 
and 3.6cr level for PSRs B1937±21 and B1957±20, respectively. Note that a fit of PSR B1937±21 
with a simple power law of the form Nq x (F/lGeV) _r gives Nq = (1.75 ± 0.29) x 10~ n cm~ 2 s" 1 
MeV' 1 and F = 3.02 ± 0.18. 


To quantify any orbital modulation of the gamma-ray flux of PSR B1957±20, we constructed a 
light curve above 0.1 GeV with 10 bins in orbital phase. We mapped each bin’s edges to a series of 
topocentric times, allowing us to correctly select photons and calculate the exposure to the source. 
We then performed a maximum likelihood analysis in which all sources except PSR B1957±20 
were held fixed at their best-fit orbit-averaged values, while the flux of the millisecond pulsar was 
allowed to vary in each bin. Its spectral shape was held fixed. We found no significant evidence 
of modulation. To place a limit on the total modulation, we first modeled the orbital modulation 
with a sinusoid. By averaging over the zero of phase, we obtained a 95% confidence upper limit on 
its peak-to-trough amplitude of 78%. However, a sinusoid is not an adequate functional form to 
characterize variability associated with only a limited range of orbital phase, so we also calculated 
the one-sided 95% confidence interval for the value of the bins with the maximum and minimum 
observed flux. To avoid bias from our choice to align the first bin with the zero of orbital phase, 
we averaged these limits over light curves shifted by 0, 1/30 and 2/30. Shifting the data by smaller 
amounts essentially had no effect on the results. The 95% limit on the maximum (minimum) flux 
so averaged is 2.1 (0.1) times the orbit-averaged flux listed in Table [2j corresponding to an upper 
limit on any excursion from the average flux of ~ 3.4 x 10~ 8 cm~ 2 s -1 . The latter value provides 
an upper limit on the emission from shock acceleration, which is consistent with the predictions of 
models of high-energy emission from colliding winds in massive stars (see e.g. Reimer et al. 2006). 


3.2.2. Gamma-ray pulsations 

With the full gamma-ray spectral models obtained from the spectral analysis of PSRs B1937±21 
and B1957±20, we were able to assign photon probabilities using the Fermi ST gtsrcprob. We 


PSR B1937+21 
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Fig. 1. - Phase- aver aged gamma-ray energy spectrum for PSR B1937+21, above 0.5 GeV. The 
black line shows the best-fit model from fitting all the gamma-ray data above 0.5 GeV with the 
simple exponentially cutoff power-law functional form given in Equation ([Tj) . Data points are 
derived from likelihood fits of individual energy bands where the pulsar is modeled with a simple 
power-law form. A 95% confidence level upper limit was calculated for any energy band in which 
the pulsar was not detected above the background with a significance of at least 2a. 


followed the prescriptions described in Kerr (2011) and calculated the weighted 17-test statistics. 
Selecting events found within 5° from the MSPs and with energies above 0.1 GeV led us to weighted 
37-test parameters of 70.4 for PSR B1937+21 and 156.1 for PSR B1957+20, corresponding to pul- 
sation significances of 7.2 and 10.9er respectively, confirming PSRs B1937+21 and B1957+20 as 
sources of pulsed gamma rays. 


Figure [3] shows the weighted light curves for the two MSPs, and as can be seen, PSR B1937+21 
exhibits two main peaks at phases ~ 0 and ~ 0.55. The error bars were derived by doing a Monte 
Carlo analysis. In each of 1000 realizations, the photon probabilities calculated with the full spectral 
model were randomly re-assigned to events in the dataset. For each, a phase histogram (a new 
weighted light curve) with the same number of bins as in Figure 3 was filled. The uncertainty for 
a given phase bin was then obtained by calculating the standard deviation of the corresponding 
phase bin contents in the shuffled light curves. If we denote Sj = w* as the probability that a 
given photon originates from the pulsar, then 6* = (1 — Wj) gives the probability that the photon 
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Fig. 2. — Same as Figure [1} for PSR B 1957+20. 


is due to background. The total background level B is obtained from the sum, weighted by bi , 
of the contribution of each photon to the weighted light curve, s,; : B = Si x 6*. The excess 
of 12.2 weighted counts above the background level observed at phase ~ 0.7 corresponds to a low 
significance of 1.3cx. On the other hand, PSR B1957+20 shows a wide emission component peaking 
at phase ~ 0.15, and a sharp ganrnra-ray peak at phase ~ 0.6. In complement to Figure [3j Figures [4] 
and [5] present phase-aligned radio and gamma-ray light curves for PSRs B1937+21 and B1957+20, 
where the gamnra-ray profiles correspond to events found within 0.8° from the pulsars and energies 
above 0.1 GeV. The absolute phasing in these light curves is such that the maxima of the first 
Fourier harmonics of the 1.4 GHz Nangay radio profiles transferred back into the time domain 
define phase 0. Under this convention, the main radio peak and interpulse for PSR B1937+21 
have their maxima at phases 4?^ ~ 0.014 and 4 >/j 2 ~ 0.536, while for PSR B 1957+20 they fall at 
phases 4?^ ~ 0.162 and 4>^ 2 ~ 0.604. The background levels above 0.1 GeV and 1 GeV shown in 
Figures [i] and [H] were obtained by calculating B = ^^(1 — w*), where w j denotes the probability 
that photon i originates from the pulsar, and N the number of photons in the considered ROI. 

We fitted the gamma-ray emission peaks shown in Figure [3] using Lorentzian functions. For 
each peak, the position 4+ and the Full Width at Half-Maximum FWHMj are listed in Table [2j 
Also listed in this Table are radio-to-gamma-ray lags <5, = 4>j — 4>/j, where 4>/j refers to the position 
of the closest radio peak. Uncertainties quoted in this Table are statistical, and are an order of 
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Pulse Phase 


Fig. 3. — Gamma-ray light curves of PSRs B1937+21 (bottom panel) and B1957+20 (top panel) 
for events recorded by the Fermi LAT within 5° from the pulsars, and with energies above 0.1 
GeV. These profiles were built by weighting each event by its probability to have been emitted by 
the pulsars, where the probabilities have been obtained from the spectral analysis of the two MSPs 
(see 3.2.1 for more details). Two rotations are shown for clarity, and the light curves have 100 bins 
per rotation. The inset shows the profile of PSR B1957+20 in the 0.57 - 0.66 phase range, with 
500 bins per rotation or ~ 3 /is per bin. 


magnitude larger than absolute phasing uncertainties due to the DM measurement (see Section 3.1 ). 
The second gamma-ray peak in PSR B1957+20’s profile is therefore found to be very sharp, with 
a width of 0.014 ± 0.007 at half maximum, or 23 ± 11 /j,s. We verified the sharpness of this peak by 
calculating the RMS of weighted phases in the 0.6 to 0.63 range. We obtained a standard deviation 
of 0.005 in phase, corresponding to a FWHM of 0.012, compatible with the FWHM measured in 
the binned analysis. This peak is the narrowest feature so far observed in a gamma-ray pulsar light 


curve. 


As can be seen from Figures [4] and [5j PSRs B1937+21 and B1957+20 show nearly aligned 


radio and gamma-ray emission components, like the Crab pulsar (Abdo et al. 2010d) and the 


millisecond pulsars PSR J0034— 0534 (Abdo et al. 2010b) and J2214+3000 (Ransom et al. 2011) 
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Fig. 4. — Multi-wavelength phase histograms of PSR B1937+21. The two bottom panels show 
radio profiles recorded at the Nangay and Westerbork radio telescopes at 1.4 and 2.3 GHz. The 
middle panel shows an X-ray light curve recorded with RXTE between 2 and 17 keV, with 100 
bins per rotation. The two top panels show 50-bin gamma-ray light curves recorded with the LAT 


above 0.1 GeV and 1 GeV. Horizontal dashed lines indicate ganuna-ray background levels (see 3.2.2 
for details on the determination of these background levels). 
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ire[4j for PSR B1957+20. The two bottom panels show radio profiles recorded 
1 Nangay radio telescopes at 0.35 and 1.4 GHz. The middle panel shows an 
. from an analysis of XMM-Newton data between 0.5 and 4.5 keV, with 25 
ie background level in this panel was extracted from a similar neighboring 
f sources. The two top panels show 100 bin gamma-ray light curves recorded 
.1 GeV and 1 GeV. 
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The radio-to-gamma-ray lags Si for the main and secondary gamma-ray peaks of PSRs B1937+21 
and B1957+20 are indeed found to be close to 0, with lcr error bars nearly as large or exceeding the 
values themselves because of the limited photon sample. In addition, for both pulsars full widths 
at half-maxima FWHMj are larger than radio-to-gamma-ray lags 6i, so that we cannot claim any 
significant separation. We conclude that for both MSPs the main and secondary emission peaks 
in gamma rays and radio seem to be produced in similar regions of the pulsars magnetosphere. 
Interestingly, while all radio peaks of PSR B1957+20 observed at 0.35 GHz have obvious gamma- 
ray counterparts, we find no evidence of a gamma-ray peak aligned with the sharp radio peak at 
phase ~ 0.95 observed at 1.4 GHz (see Figure [5]), suggesting strong spectral dependence of the 
emission regions and correlation between the gamma-ray and low-frequency radio emission. 


With the firm detection of pulsed gamma-ray emission from PSR B1957+20, we can now 
identify the millisecond pulsar as the source powering the Fermi LAT First Year Catalog source 
1FGL J1959. 6+2047, located 0.5' away (Abdo et al.| 2010c). The energy flux above 0.1 GeV 
measured for the pulsar in this analysis (see Table [2]) is in agreement with that of the 1FGL 
source, of ~ 2.0 X 10~ n erg cm~ 2 s . In the case of PSR B1937+21, two gamma-ray sources 
with no known counterparts are found within a few tens of arc minutes: 1FGL J1938.2+2125c and 
J1940.1+2209c, with best-fit positions located 21.7’ and 35.6’ away, respectively. Although the 
energy flux for PSR B1937+21 above 0.1 GeV is formally consistent with that of the latter 1FGL 
source of ~ 3.4 x 10 -11 erg cm -2 s , it was still detected as a separate ganmra-ray source from 
the pulsar B1937+21, while the former was no longer significant. This, in addition to the relatively 
large angular separations and the fact that the 1FGL sources are flagged as being possibly spurious 
or confused with the diffuse emission, indicate that the MSP probably is not formally associated 
with either of the 1FGL sources. Instead, they may result from the source confusion due to the 
presence of a weak gamma-ray pulsar in a region of intense background. 


3.3. X-ray analysis 

3.3.1. PSR B 1937 +21 


To study the relative alignment of X-ray and gamma-ray emission components, we have re- 
analyzed the data presented in Cusumano et al. (2003). Observations of PSR B1937+21 were 
made using the Proportional Counter Array (PGA; Jahoda et al. 2006) on board the Rossi X-ray 
Timing Explorer ( RXTE ). The PCA is an array of five collimated Xenon/methane multi-anode 
proportional counter units (PCUs) operating in the 2-60 keV range, with a total effective area of 
approximately 6500 cm 2 and a field of view of ~ 1° FWHM. Data were collected in “GoodXenon” 
mode, which records the arrival time (with 1 //s resolution) and energy (256 channel resolution) 
of every un-rejected event. Twenty-one observations taken between 2002 February 21 and 2002 
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February 28 were downloaded from the HEASARC archive]^] and standard selection criteria were 
applied. Photon arrival times were converted to barycentric dynamical time (TDB) at the solar 
system barycenter using the JPL DE200 solar system ephemeris with the FITS tool faxbary. In 
order to achieve the best possible absolute timing accuracy available with RXTE of ~ 10 /is (e.g. 
Rots et al. 1998, including the uncertainty introduced from conversion between TT and TDB), the 


fine clock correction was applied to each evenlj^] In order to be consistent with the previous analysis 
of these data described in Cusumano et al. (2003), we selected photons from all three xenon layers 
in the energy range 2-17 keV. This analysis produced a total of 4.46 X 10 6 photons. Using the 
parameters and absolute phase information provided in Cusumano et al. (2003), and correcting an 
error in the value for the derivative of the DM which was incorrectly given as 2.1 pc cm -3 yr” 1 
instead of 2.1 x 10” 3 pc cm” 3 yr” 1 (M. Kramer, private communication), we folded all of the X-ray 
photons and created a single pulse profile, shown in Figure [4] 

Our results are consistent with those of Cusumano et al. (2003): a fit of the two X-ray emission 
peaks places the first component at 4 > x 1 = 0.0512 ± 0.0003 with a FWHM of 0.018 ± 0.001 and the 
second one at (f>x 2 = 0.5752 ± 0.0004 with a FWHM of 0.040 ±0.012, in phase units. Error bars are 
statistical, and do not account for the timing accuracy of RXTE of ~ 10 /rs, or 0.006 in phase units. 
We therefore confirm that the X-ray peaks of PSR B1937±21 coincide with the giant radio pulse 
emission and not with the regular radio emission. Additionally, from the positions and widths of the 
gamma-ray peaks listed in Table [2j it seems that X-ray and gamma-ray emissions are misaligned, 
indicating that they are produced in slightly different regions of the magnetosphere. Note however 
that the DM variations reported in Cusumano et al. (2003) were not observed in the radio timing 
analyses made in support of Fermi observations (see Section 3.1). In order to exclude any relative 
phasing issues induced by wrong DM estimates or instrumental issues, we have undertaken the 
analysis of other archival data supported by multi-frequency radio timing observations spanning 
over several years. This work will be presented in a future paper. 


3.3.2. PSR B 1957 +20 


No detection of X-ray pulsations have been reported for PSR B1957±20 thus far (see Huang 


&: Becker] 2007. for a recent analysis). However, it must be noted that folding the ganmra-ray data 


for PSR B1957±20 with the pulsar ephemeris used in the latter paper does not reveal gamma-ray 
pulsations. While this shows that the ATNF Catalogue timing solution used in Huang & Becker 


(2007) is invalid for the Fermi data taken after August 2008 it does not prove that the ephemeris 
was invalid for the XMM-Newton data taken in 2004. This motivated a re-analysis of the same 


X-ray data, using the ephemeris presented in Section 3.1 


7 http:/ /heasarc. gsfc.nasa.gov/docs/archive.html 
8 ftp: / /heasarc. gsfc/xte/calib_data/clock/tdc.dat 


16 - 


PSR B1957+20 was observed with XMM-Newton on 2004 October 31 and 2004 November 1. 
The observations with the three European Photon Imaging Cameras (EPIC) MOS1, MOS2 (imaging 
mode) and pn (timing mode) spanned approximately 30 ks. These instruments were operated in 
full-frame mode with a thin filter. We used Version 10.0 of the XMM-Newton Science Analysis 
Software]^] to reduce and analyze the X-ray data. The EPIC Observation Data Files (ODFs) were 
reduced using the emproc/epproc scripts (for MOS and pn respectively) along with the most recent 
Calibration Files (CCFs) at the time of the reduction. We applied standard filtering procedures 
( Web b~et al.|| 2004 ) , which includes the removal of strong background periods caused by soft photon 
flares, and considered the data between 0.3 and 10 keV for the spectral analysis. We restricted 
the timing analysis to the 0.5 - 4.5 keV band as this was found to have the best signal-to-noise. 
We removed the events contained in energy ranges affected by internal background caused by the 
X-ray fluorescence of satellite material exposed to cosmic rays 


10 


We extracted MOS spectra from circular regions of radii 45 w centered on the radio position 
of the pulsar, which allowed us to collect about 90% of all detected source counts and optimize 
the signal to noise ratio. In the pn timing mode, the central CCD is read out continuously at 
high speed, and the collected events are condensed into a one dimensional pixel array. The source 
event file is therefore extracted by selecting a rectangular region around the expected position of 
the pulsar. We extracted MOS and pn background photons from neighboring regions free of X-ray 
sources, and generated instrumental response files with the RMFGEN/ ARFGEN tasks. 

We fitted the combined MOS1/2 and pn spectra with X 'spec Version 12.6. <E3 Using a simple 


power-law model yielded spectral results that are consistent with those of Huang & Becker (2007), 


with a photon index T = 2.37+0 29 (x 2 = 0.70 for 27 d.o.f., errors are 90%). The fitted column 


absorption along the line-of-sight Nh is 15.7_ 9 ' 0 x 10 2U cm , and is consistent with the total 
Galactic HI column density in the direction of the pulsar given by the HEASARC Tool Vh/F^l 
(~ 30 X 10 20 cm -2 ). The unabsorbed X-ray flux in the band 0.2 - 10 keV is F = (9.7 ± 0.8) x 10 -14 
erg cm -2 s -1 . Using the NE2001 distance of 2.5 kpc, the X-ray luminosity is found to be Lx = 
(7.2 ± 0.6) x 10 31 erg s -1 . Assuming the emission is produced by the pulsar, we estimate that 
the efficiency of the conversion of spin-down energy loss E into X-ray emission is rjx = Lx/E ~ 
9.7 x 10 -4 . These values are higher than the ones presented by Huang & Becker (2007), who 


considered a distance to the pulsar of 1.5 kpc. However, with the lower limit on the distance of 


d> 2 kpc placed by van Kerkwijk et al. (2011), the X-ray luminosity and efficiency values measured 


in this analysis are favored. 


The relative timing accuracy of XMM-Newton has been proven to be better than 10 8 , while 


9 http:/ /xmm.esac.esa.int/sas/ 

10 http:/ /xmm2.esac.esa.int/docs/documents/C AL-TN-0018.pdf 
1 1 http : / /heasarc .gsfc.nasa. gov / xanadu /xspec/ 

12 http:/ /heasarc. gsfc.nasa.gov/cgi-bin/Tools/w3nh/w3nh. pi 
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the absolute timing accuracy is ~ 53 /as (Martin-Carrillo et al. 2011), allowing to search for pul- 
sations in the X-ray data using the pn source event file extracted as described above. The event 
times were converted to TDB using the task barycen Version 1.18 and the JPL DE405 solar sys- 
tem ephemeris, and phase-folded using Tempo2 and the ephemeris described in Section 3.1 This 


timing solution however does not formally cover the X-ray data considered here, taken four years 
before the beginning of the radio timing dataset. Millisecond pulsars are generally stable rotators, 


therefore it is reasonable to extrapolate the ephemeris to our observations in a similar way to Bog- 


danov et al. (2010). Arzoumanian et al. (1994) measured a DM for PSR B1957+20 of 29.1168 db 


0.0007 pc cm 3 at 48196 MJD. We can therefore safely assume that the DM value at the time 
of the XMM-Newton observations was comprised between the latter value of Arzoumanian et ah 


(1994) and our DM measurement of 29.12644 ± 0.00029 pc cm' 3 at 54550 MJD. The difference 
in DM between 48196 and 54550 MJD of ~ 9.6 x 10 -3 pc cm' 3 introduces an uncertainty on the 
relative phasing of XMM-Newton data and 1.4 GHz radio data of ~ 20 [is. This relatively small 
value indicates that the uncertainty on the DM at the time of the X-ray observations should not 
affect our timing analysis of the XMM-Newton data significantly. Nevertheless, our analysis could 
be affected by the instability of the pulsar’s rotational frequency or orbital movement. 


The X-ray phase histogram is shown in Figure [5] The pulse profile gathers 854 events taken 
between 0.5 and 4.5 keV, of which ~ 33% are estimated to come from the pulsar. The H - test 
parameter for this set of events is 24, which corresponds to a pulsation significance of ~ 4a. 
Assuming the ephemeris was accurate at the time of the XMM observations, it hence seems that 
PSR B1957+20 emits pulsed X-rays. The background level shown in Figure [5] was calculated by 
extracting a similar neighboring region free from X-ray sources, as mentioned above. We fitted 
the X-ray peak with a Gaussian and found that it occurs at (f>x = 0.86 ± 0.03, with a FWHM 
of 0.32 ± 0.03. It is therefore found to be offset from the aligned radio and gamma-ray peaks 
by at least 150 fis. From the absolute timing uncertainty of ~ 73 //s ~ 0.05 in phase induced by 
XMM-Newton ' s intrinsic absolute timing capability and the uncertainty on the DM, it seems that 
the X-ray peak is separated from the phase-aligned radio and the gamma-ray emission peaks. We 
cannot exclude however drifting induced by non-optimal rotational and binary parameters, causing 
the pulse phases to be offset from the actual values. A longer X-ray observation of the MSP, 
analyzed with a contemporaneous timing solution would be preferable for confirming the phase and 
shape of the peaks. 


Discussion 


4.1. Light curve modeling 


The first eight MSPs discovered using the Fermi LAT (Abdo et al. 2009a) exhibited non-zero 
lags between their respective radio and gamma-ray light curves, and have been modeled using 
either standard two-pole caustic (TPC) and outer gap (OG) geometries, or pair-starved polar cap 
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(PSPC) models (Venter et al.||2009), or alternatively, an annular gap model (Du et alT||2010 ) . Light 


curve modeling using a force- free magnetospheric geometry may present a further possibility (Bai 


&; SpitkovskypOlO Contopoulos fc Kalapotharakos||2010 ) . The common feature of all these models 
is that the gamma-ray emission originates in the outer magnetosphere and not near the polar caps. 


The discovery of radio and gamma-ray light curve peaks in close alignment for J0034— 0534 (Abdo 


et al. 2010b) yielded the first example of an MSP with such phase-aligned peaks, previously only 


seen 


in the Crab pulsar. The newly-discovered PSR J2214+3000 (Ransom et al. 2011) also shows 
aligned radio and gamma-ray emission peaks. As can be seen in Figures [4] and [5] showing multi- 
wavelength light curves for PSRs B1937+21 and B1957+20, we now have a class of MSPs with 
phase-aligned low and high-energy pulsations, in addition to PSRs J0034— 0534 and J2214+3000. 

Light curve modeling efforts for PSRs B1937+21 and B1957+20 are guided by the following 
key observed features: (i) gamma-ray and radio peak phase alignment within the statistical uncer- 
tainties; (ii) gamma-ray profiles that have multiple sharp peaks and possible increased complexity 
if additional low-level features become more significant with accumulated photon statistics. 


As described in Abdo et al. (2010b), the alignment of gamma-ray and radio peaks suggests 


co- location of gamma-ray and radio emission regions in the pulsar magnetosphere. We obtain 
reasonable fits for both gamma-ray and radio light curves in the context of “altitude-limited” TPC 
(alTPC) and OG (alOG) models (Figures [6] and [7j . These models assume that the radio emission is 
emitted within a range of altitudes relative to the light cylinder (at radius Rlc = cP/2ir), along the 
last open field lines, with the same geometry as the gamma-ray emission except that the emission 
is limited to a smaller range of altitudes. Presently, high-altitude emission seems to be preferred 
to low-altitude emission, although the latter cannot be ruled out (Venter et al. 2011). Thus, 


both the gamma-ray and radio peaks are caustics, resulting from phase-bunching of photons due 


to relativistic effects associated with high corotation velocities in the outer magnetosphere (Dyks 


et al.|2004 Romani fc Yadigaroglu|1995 ). For most viewing geometries, caustic radio emission leads 
to large amounts of depolarization due to mixing of emission from different altitudes; however, 
orthogonal configurations with viewing angles near 90° result in less depolarization and thus a 
careful study of the expected polarization properties in these models is warranted. 

In order to statistically pick the best-fit parameters for the alTPC and alOG models for the 
two MSPs considered here we have developed a Markov chain Monte Carlo (MCMC) maximum 
likelihood procedure (Johnson et al.||2011). An MCMC involves taking random steps in parameter 


space and accepting a step based on the likelihood ratio with respect to the previous step (Hastings 


1970). The gamma-ray light curves are fitted using Poisson likelihood while the radio light curves 
are fitted using a x 2 statistic and the two values are combined. For a given parameter state 
the likelihood value is calculated by independently optimizing the radio and gamma-ray model 
normalizations using the scipy python modulcp*| and the scipy . optimize. fminJ-fbgs-b multi- variate 
bound optimizer (Zhu et al. 1997). The likelihood surfaces can be very multi-modal which can lead 


See http://docs.scipy.org/doc/ for documentation. 
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to poor mixing of the chain and slow convergence. Therefore, we have implemented small-world 
chain steps (Guan et al. 2006) and simulated annealing (Marinari & Parisi 1992) to speed up the 
convergence and ensure that the MCMC fully explores the parameter space and does not get stuck 
in a local maximum. We verify that our chains have converged using the criteria proposed by 


Gelrnan & Rubin (1992). 


In order to balance the gamma-ray and radio contributions to the likelihood, we have chosen 
to use an uncertainty for the radio intensity which is equal to the average, relative gamma-ray 
uncertainty in the on-peak region times the maximum radio value. The best-fit results can be 
strongly affected by the uncertainty which is chosen for the radio profile; in particular, a smaller 
uncertainty will decrease the overall likelihood and can, in some cases, lead to a different best-fit 
geometry which favors the radio light curve more strongly. When varying the radio uncertainty by a 
factor of 2, the best-fit a and ( values of PSR J1939+2134 were found to vary by < 7°. The best-fit 
geometry of PSR J1959+2048 was found to be more sensitive to changes in the radio uncertainty 
with either the best-fit aor( value changing by ~ 35° while the other parameter changed by < 15°. 


Our simulations have a resolution of 1° in both a and £, 0.05 in gap width (w, normalized 
to the polar cap angle 9 pc ~ yj R ns/ Rlc), and 0.1 Rlc for the minimum and maximum emission 
altitudes. The MCMC explores viewing geometries for a from 1° to 90° and for ( from 0° to 
89°, as shown in Figure [8} We expect that going beyond 90° should produce the same profiles, 
simply shifted in phase by 180°. Additionally, for the gamma-ray models we restrict the maximum 
emission altitudes to be > 0.7 Rlc, as high-altitude emission near the light cylinder is important 
for producing the correct pulse shapes. In addition to predicting the best- fit model parameters the 
simulations also provide numerical estimates of the beaming correction factor (/q) as described 
in Watters et al. (2009) and Venter et al. (2009). For a given a and set of emission parameters 


(i.e. , gap width and emission altitudes), the /q factor is calculated by collecting the emission at a 
specified ( and comparing it to the total emission for all viewing angles. 


The radio and gamma-ray profiles of PSRs B1937+21 and B1957+20 were reproduced using 
the alTPC and alOG geometries with the best-fit parameters listed in Table [3j the uncertainties 
on the model parameters were derived from a likelihood profile scan. In this Table, Rns denotes 
the neutron star radius, and Rncs is the field-line-dependent altitude of the null-charge surface, 
defined by 12 • B = 0. Best-fit models with infinitely thin gap widths (w = 0) do not represent 
the truth as a zero-width gap is unphysical. However, they indicate that the best gap width is 
somewhere between 0 and 0.05 and the best-fit value of 0 is chosen only as a result of the resolution 
of our simulations. For PSR B1937+21, the fits yield large a and ( values, which reinforces the idea 
that this MSP may be a nearly orthogonal rotator, as is suggested by the fact that the observed 
radio iuterpulse lags the main radio peak by approximately half a rotation. Recent observations 
using the Parkes telescope indicate that PSR B1937+21 may have a large duty cycle of ~ 80% 
(Yan et al. 2011). The two main peaks separated by 172° remain the main features of the profile, 
while off-peak features at ~ 0.5% of the peak intensity are seen. While the current radio models 
reproduce the two main peaks very well, it may be worthwhile to model the extended low-level 
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off-peak emission as well in future, as these features may result from emission regions below the 
null charge surface. If this is true, emission from both poles is implied, favoring a TPC geometry. 


The alTPC and alOG models also yield similar values of the a and ( angles for PSR B1957+20, 
preferring moderate values of a and ( near 90° . The spin-up of millisecond pulsars by mass transfer 


from a companion star naturally leads to aligned orbital and rotation axes (Thorsett Sz Stinebring 


1990). Since the pulsar shows eclipses (Fruchter et al. 1988), the line-of-sight is nearly parallel to 


the plane of the orbit, which indicates that a predicted observer angle £ close to 90° is preferred. 
However, neither model is able to successfully reproduce the narrow gamma-ray peaks well and, 
as can be seen in Table [3j a is not well constrained, particularly for the alOG model. With more 
statistics the gamma-ray peaks will be more important to the likelihood and the fits may improve. 


4.2. Constraints on orientation angles from radio data 


In contrast with most other MSPs, the radio polarization profile of PSR B1937+21 is simple, 
with a distinct flat polarization angle (PA) swing (see e.g. Stairs et al. 1999). This suggests that the 
orientation angles of the pulsar can be inferred by fitting the polarization data with the rotating 
vector model (RVM, Radhakrishnan & Cooke 1969), accounting for the aberration effect bending 
the emission beam in the co-rotational frame of the pulsar forward with respect to the rotating 
frame, and resulting in a delay of the PA swing respective to the intensity profile by Ar/R^c radians, 
where r is the emission radius (Blaskiewicz et al. 1991). We therefore attempted to model the radio 
emission geometry, using the 0.6 and 1.4 GHz polarization data from Stairs et al. (1999), obtained 
from the EPN database 14 The degree of linear polarization of the main peak is fairly high (> 50%) 
while that of the secondary peak is 10%. A PA jump due to orthogonal polarization modes (OPM) 
is identified in the main peak at 1.4 GHz, and in the secondary peak at 0.6 GHz. The continuity 
of the OPM jumps and the fact that they are not exactly 90° apart may imply that the PA swings 
are mildly affected by scattering in the interstellar medium (Karastergiou 2009). 


The PA profiles at 0.6 and 1.4 GHz corrected for the OPM jumps are very similar. In addition, 
RVM fits of the two PA profiles give consistent results. We therefore merged the two profiles into 
a single PA swing. Figure [8] shows the results of the RVM fit to the merged PA swing. The best fit 
of the data is obtained for a and /3 angles of 89° and —3° respectively, where f3 = C ~ a is the angle 
between the magnetic axis and the line of sight. Also shown in the plot are the pulsar orientation 
angles obtained from the modeling of the radio and gamma-ray profiles under the alTPC and alOG 
geometries. As can be seen from Figure [8j the best-fit configuration is found in the vicinity of 
the angles obtained from the modeling. Both the radio and gamma-ray light curve modeling and 
the analysis of the polarization data therefore support values of the a and C angles close to 90°, 
corresponding to an orthogonal configuration. 


14 http:/ /www. mpifr-bonn.mpg.de/div/pulsar/data/browser.html 
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Fig. 6.— Top: gamma-ray data and modeled light curves for PSR B1937+21 with 30 bins per 
rotation. Bottom: Nangay 1.4 GHz radio profile and modeled light curves. See Table [3] for the 
best-fit parameters. 



Fig. 7. — Same as Figure [6j for PSR B1957+20. The bottom panel shows the Westerbork 0.35 
GHz radio profile. 
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Fig. 8. — Plot of the impact angle /3 = ( — a, vs. the magnetic inclination angle a for 
PSR B1937+21. The y 2 map of the rotating vector model (RVM) fit of radio polarization data for 
PSR B1937+21 is shown in gray scale, with contour levels at 1, 2 and 3 a in black. Magenta and 
green shaded regions indicate the 3er contour levels of the simultaneous radio and ganrma-ray light 
curve modeling under the alTPC and alOG geometries, while the solid lines show the 1 a contour 
levels. The symbols indicate the best pulsar orientation angles obtained from the RVM, alTPC 
and alOG fits. 
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In addition to inferring pulsar geometric angles, the RVM fit also allowed us to estimate the 
altitude of the radio emission, and compare it to the modeled values listed in Table [8} We find that 
the magnetic axis is located at ~ 150° after the secondary peak (~ 0.42 in phase), corresponding to 
a radio emission altitude of 0.65 Rlc- This value is consistent with the emission altitude found from 
the modeling under both alOG and alTPC geometries, confirming that the radio emission seems 
to be produced at high altitudes above the neutron star surface. We note that the aberration 


treatment presented by Blaskiewicz et al. ( 1991 ) and later by Hibschman & Arons (2001 ) and Dyks 


(2008) is only a first-order approximation that may need correcting, in particular at high altitudes. 


However, the consistency of the results obtained from the RVM fitting and the light curve modeling 
gives us some confidence in suggesting that both radio and gamma-ray emission originate from the 
outer magnetosphere, which for PSR B1937+21 is compact anyway. 

On the other hand, radio polarization data for PSR B1957+20 could not be used to constrain 
the pulsar orientation angles. This pulsar displays less than 2% linear polarization though some 


evidence exists for sign-changing, circular polarization through both peaks (Thorsett Sz Stinebring 


1990). Observations well away from the eclipsing phase have confirmed that the lack of polarization 


is not due to interactions with the companion. This lack of linear polarization is consistent with 
radio emission of a caustic nature. Note that similar behavior is observed for PSR J0034— 0534 
(|Stairs et al. 1999) which has also been fit with the alTPC and alOG models (Abdo et al. 2010b 


Venter et al. 2011). Future modeling efforts should be able to reproduce the weak polarization of 


caustic peaks. 


4.3. Gamma-ray luminosities 

Knowing the pulsar distance d and the gamma-ray energy flux G measured using the LAT above 
0.1 GeV, one can derive the gamma-ray luminosities L 7 = Anf^Gd 2 and efficiencies of conversion of 
spin-down energy into gamma-ray emission r] = L^/E of PSRs B1937+21 and B1957+20. In these 
expressions, /q is the correction factor depending on the viewing geometry discussed in Section 
|4.1| Table [2] lists values of L 7 and r] under the assumption that fn = 1 . The geometrical correction 
factors /n obtained from the modeling of radio and gamma-ray light curves with alTPC and alOG 
geometries are listed in Table [3j as well as the corresponding corrected gamma-ray luminosity and 
efficiency values. 

With gamma-ray luminosity estimates on the order of ~ 2.5 x 10 35 erg s“\ PSR B1937+21 
does not stand out from gamma-ray pulsars with comparable E values (see e.g. Abdo et al. |]2010e) . 
Nevertheless, the gamma-ray luminosity above 0.1 GeV inferred from using the NE2001 distance 
of 3.6 ±1.4 kpc is L 7 / fn = (5.63 ± 3.95 ± 3.53) x 10 34 erg s _1 (where the first uncertainty quoted is 
statistical and the second is systematic) . This value is consistent with those of pulsars with similar 
E values within uncertainties, and therefore neither parallax distance estimates nor distances based 
on Galactic electron density models are favored with the current gamma-ray analysis and knowledge 
of the relationship between ± 7 and E. 
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Similarly, the gamma-ray luminosities measured for PSR B1957+20 are consistent with those 
of pulsars with comparable E values to within uncertainties, despite the fact that the /q correction 
factors obtained from the modeling are very different between the two models. We therefore cannot 
constrain the distance of 2.5 kpc inferred from the NE2001 model of Galactic electron density. 
Nevertheless the low gamma-ray efficiency inferred from the alOG model is very similar to those 
of the bulk of MSPs (Abdo et al. 2009a). This fact, in addition to the inferred £ being close to 
90° as suggested from the occurrence of radio eclipses and the slightly better likelihood parameter 
suggests that the altitude-limited OG model is preferred over the TPC model for PSR B1957+20. 


5. Summary and conclusions 


We have reported the discovery of phase-aligned radio and gamma-ray light curve peaks for 
the millisecond pulsars PSRs B1937+21 (the “first” MSP) and B1957+20 (the first black widow 
system). This adds two new members to the class of MSPs exhibiting the phenomenon of such 
phase-aligned peaks. 

The fact that we find reasonable (but possibly non-unique) radio and gamma-ray light curve 
fits implies that the geometric caustic models still provide an adequate description for this new 


class of MSPs. As noted before (Venter et al. 2009), the sharp peaks, coupled with the caustic 


fits, imply copious pair production in the MSP magnetospheres, in contrast to earlier expectations. 
Cascades of electron-positron secondaries are needed to set up and sustain the TPC / OG geometries 
which reproduce the salient features of the light curves. PSR B1937+21 is very near the death 
line for screening by pairs, but PSR B1957+20 lies well below the line and is not expected to 
produce enough pairs for screening in conventional polar cap models assuming dipolar magnetic 
fields (Harding et al. 2005). In addition, these geometric models provide a framework to constrain 
the emission altitudes of the gamma-ray and radio photons, as well as the beaming and inclination- 
observer geometry of the MSPs. Lastly, observations of the radio emission at different frequencies, 
as well as energy-dependent light curve modeling, may provide the opportunity to learn more about 
the radius-to-frequency mapping of the radio, its connection with the gamma-ray radiation, and 
ultimately the mechanism for the generation of the radio emission. 

Our analysis of the RXTE X-ray data for PSR B1937+21 yielded results that are consistent 


with those of Cusumano et al. ( 2003 ) . The X-ray light curve consists of two peaks lagging the regular 
radio emission by a small amount, but in close alignment with the giant radio pulse emission. X- 
ray pulses also are not formally aligned with the gamma-ray emission seen with the LAT. Detailed 
modeling of the X-ray emission geometry for this MSP would help understand the misalignment 
with the regular radio emission and the gamma-ray emission. Nevertheless, the sharpness of the X- 
ray peaks and their proximity to the outer magnetospheric radio and gamma-ray emissions suggest 
that the X-ray emission from PSR B1937+21 also takes place at high altitude in its magnetosphere, 
which is supported by the non-thermal nature of the emission ( Cusumano et al.|2003 Nicastro et al 
20041). 


- 25 - 


On the other hand, we have found evidence for X-ray pulsed emission from PSR B1957+20 
using XMM-Newton data, for the first time. The current timing analysis suffers from potential 
phase drifting due to the fact that our radio timing solution does not cover the X-ray data epoch. 
Nevertheless, the relatively large spin-down luminosity E of this pulsar of 7.5 x 10 34 erg s _1 , 
comparable to that of other X-ray emitting MSPs, makes it a good candidate for pulsed X-ray 
emission. Besides, this spin-down luminosity is characteristic of non-thermally emitting MSPs, 
which is supported by the spectral analysis of the XMM-Newton data. If PSR B1957+20 does emit 
pulsed X-rays with a non-thermal spectrum, then the X-ray emission from this MSP is expected 
to align with the radio emission, as is the case for PSRs B1937+21 and J0218+4232 (Kuiper et al. 


2002). Is is therefore important to observe this pulsar in X-rays again, with a radio timing solution 


covering the X-ray data and ensuring the validity of the absolute phasing. 
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Table 1: Properties of PSRs B1937+21 and B1957+20. Values in parentheses indicate the ler 
uncertainties on the last digit quoted. For PSR B1937+21, P and P and ht values are taken from 


Cognard et al. 

(1995 

), while values for PSR B1957+20 are taken from 

Arzoumanian et al. 

(1994) 


These value are given at epochs MJD 47899.5 and 48196 in TDB units, respectively. Distances d 
are derived from timing parallax measurements from Verbiest et al. (2009) for PSR B1937+21, and 
from the NE2001 model of Galactic electron density (Cordes & Lazio 2002) for PSR B1957+20. 
The total apparent proper motion for PSR B1937+21 is small, making the Shklovskii contribution 


to the period derivative value (Shklovskii 1970) almost negligible. The last three parameters were 


calculated using the intrinsic spin-down rate P CO m corrected for the Shklovskii effect. 


Parameter 

PSR B1937+21 

PSR B1957+20 

Galactic longitude, l (deg) 

57.51 

59.20 

Galactic latitude, b (deg) 

-0.29 

-4.70 

Pulsar period, P (ms) 

1.557806472448817(3) 

1.60740168480632(3) 

Apparent period derivative, P (10~ 21 ) 

105.1212(2) 

16.8515(9) 

Transverse proper motion fj, x (mas yr _1 ) 

0.80(2) 

30.4(6) 

Distance d (kpc) 

7.7 ± 3.8 

2.5 ± 1.0 

Corrected period derivative, P CO rr (10~ 21 ) 

105.10 ± 0.01 

7.85 ± 3.61 

Spin-down luminosity, E (10 34 erg/s) 

109.76 ± 0.01 

7.48 ± 3.43 

Surface magnetic field, B sur f (10 s G) 

4.0946 ± 0.0003 

1.12 ± 0.52 

Light cylinder magnetic field, B L q (10 5 G) 

9.8472 ± 0.0006 

2.49 ± 0.81 


Table 2: Light curve and spectral parameters of PSRs B1937+21 and B1957+20 in gamma rays. 
Details on the measurement of these parameters are given in Section 3.2 Peak positions full 
widths at half-maxima FWHMj and radio-to-gamma-ray lags 5i are given in phase units, between 
0 and 1. 


Parameter 

PSR B1937+21 

PSR B1957±20 

First peak position, 4>i 

0.004 ± 0.009 

0.146 ±0.026 

First peak full width at half maximum, FWHMi 

0.030 ±0.029 

0.137 ±0.074 

First peak radio-to-gamma-ray lag, <5i 

-0.010 ±0.009 

-0.016 ±0.026 

Second peak position, $2 

0.543 ±0.013 

0.616 ±0.002 

Second peak full width at half maximum, FWHM 2 

0.041 ±0.041 

0.014 ±0.007 

Second peak radio-to-gamma-ray lag, 82 

0.006 ±0.013 

0.012 ±0.002 


Photon index, T 



1.43 

± 

0.87 

± 

0.40 

1.33 

± 

0.57 

± 

0.09 

Cutoff energy, E c (GeV) 



1.15 

± 

0.74 

± 

0.43 

1.30 

± 

0.56 

± 

0.13 

Photon flux, F (> 0.5 GeV) (10 -i 

3 cm 2 s x ) 


1.22 

± 

0.23 

± 

0.05 

0.77 

± 

0.09 

± 

0.01 

Energy flux, G (> 0.5 GeV) (10 -1 

1 1 —2 
erg cm s 


1.98 

± 

0.32 

± 

0.04 

1.34 

± 

0.15 

± 

0.01 

Extrapolated photon flux, F (> 0 

.1 GeV) (10~ 8 

cm 2 s 4 ) 

5.97 

± 

4.89 

± 

3.58 

3.09 

± 

1.62 

± 

0.43 

Extrapolated energy flux, G (> 0. 

1 GeV) (10- 11 

erg cm -2 s _1 ) 

3.63 

± 

1.58 

± 

1.09 

2.17 

± 

0.54 

± 

0.11 

Luminosity, ± 7 / /n (> 0.1 GeV) 

(10 34 erg s _1 ) 


25.8 

± 

21.2 

± 

19.6 

1.62 

± 

1.00 

± 

0.92 

Efficiency, r) / fa (> 0.1 GeV) 



0.23 

± 

0.19 

± 

0.18 

0.22 

± 

0.17 

± 

0.16 
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Table 3: Best-fit parameters obtained from the modeling of radio and gamma-ray light curves of 
PSRs B1937+21 and B1957+20, as discussed in Section 4.1 , for each emission geometry, altitude- 
limited Two Pole Caustic (alTPC) or altitude-limited Outer Gap (alOG). Altitudes are expressed 
relative to the light cylinder radius. See the text for more details on the different parameters. 



PSR B1937+21 

PSR B1957±20 

Parameter 

alTPC 

alOG 

alTPC 

alOG 

Magnetic inclination angle, a (°) 

75t? 

84t^ 

47+ 5 „ 

31 I 3 9 

Observer angle, £ (°) 

80 ± 3 

84^3 

85±1. 

89±1 

Gamma-ray emission gap width, w 7 

0.10 ±0.05 

0.05 ±0.05 

0.05 ±0.05 

0.05 ± 0.05 

Radio emission gap width, wr 

0.00 ±0.05 

0.00 ±0.05 

0.05 ±0.05 

0.10 ±0.05 

Gamma-ray emission altitudes 

[ Rns '■ 1 ± 0.2] 

[«jvcs; 1±8:?] 

[fl ws ; 1 . 2 lg;i] 

[fijvcs: l-ltoi] 

Radio emission altitudes 

[o.7±g;|; 0.9±g;?l 

[0.6 ±0.1; 0.9 ± 0.1] 

1 0.8 ± 0.1; l.O^o'iJ 

[0.7 ±0.1:0.91^1 

Geometrical correction factor, /q 

-l n +0.08 

±-U — 0.03 

n os - ! -9 - 93 

U.9o_o 02 

U.OD-o 02 

n 00+0.O6 
U ‘°^-0.12 

Likelihood parameter, — ln(L) 

126.3 

130.9 

123.7 

128.3 

Corrected L 7 (> 0.1 GeV) (10 34 erg s _1 ) 
Corrected 77 (> 0.1 GeV) 

Qr q + ^ 1 . 8 +] 2 U .]2 

zo.o _ 24 8-23.4 

0 qo~ 1 - 0 . 20 -I- 0.18 

u ‘ zc) -0.23-0.21 

q c 9 + 2 i.u+iy .4 
zo,z _27. 1-26.0 
0 q 0 — (— 0.19 — 1 ~ 0.18 
U - ZO -0.25-0.24 

Q iy-i 1 .&'Z- pl. 7 U 
‘ -*-- 2 . 22 - 2.12 
n q«~I-0.29-|-0.28 
U,OU -0.34-0.33 

n 0*7+1.64-1-1.63 

u - 0< -0.35-0.33 

r\ r\r + 0 . 22 + 0.22 

u,uo — 0.05— 0.05 


